///////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                       //
//  ------------------------------- Ebrahim Foulaadvand, 11 may 2012 ----------------------------------  //
//                                                                                                       //
//    The programme "DoubleNuclearDecay" evaluates the time dependence of the number of undecayed        //
//    radioactive nuclei of types A an B. Comparison is done with analytical results.                    //           //
//                                                                                                       //
//                                                                                                       //
//                                                                                                       //
///////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>


using namespace std;


double  tau=0.01,rA=0.2,rB=0.3; // tau=time step, rA=nucleus A decay rate, rB=nucleus B decay rate.

int T=2500,N0A=1000,N0B=200,n; // T=number of time steps.

// N0A=intial number of radioactive nucleus A, N0B=intial number of radioactive nucleus B

main()

{

ofstream  fileNA ("NA-t rA=0.2 rB=0.3.plt"); //output file for the computed time dependence of NA
ofstream  fileNAanal ("NAanal-t rA=0.2 rB=0.3.plt"); //output file for the analytic time dependence of NA
ofstream  fileNB ("NB-t rA=0.2 rB=0.3.plt"); //output file for the computed time dependence of NB
ofstream  fileNBanal ("NBanal-t rA=0.2 rB=0.3.plt"); //output file for the analytic time dependence of NB


vector<double> NA(T+1,0),NAanal(T+1,0),NB(T+1,0),NBanal(T+1,0);

// Arrays NA and NB store the computed number of undecayed nuclei A and B at time step n.
// Arrays NAanal and NBanal store the analytic number of undecayed nuclei A and B at time step n.


NA[0]=N0A,NB[0]=N0B;  // ----------------------- initial conditions --------------------


// ----------------------------- Euler algorithm   ----------------------------------------

for(int n=0;n<T; n++){

NA[n+1]=NA[n]*(1-rA*tau);

NB[n+1]=NB[n]*(1-rB*tau) + rA*tau*NA[n];


fileNA<<n*tau<<" "<<NA[n]<<endl;

fileNB<<n*tau<<" "<<NB[n]<<endl;

fileNAanal<<n*tau<<" "<<N0A*exp(-rA*n*tau)<<endl;

fileNBanal<<n*tau<<" "<<((rA/(rB-rA))*N0A*(exp(-rA*n*tau)-exp(-rB*n*tau)) + N0B*exp(-rB*n*tau))<<endl;

}

cout<<"end."<<endl;getch();

return 0;

}